### Political fragmentation over time
### 1 March 2023


# Data preparation --------------------------------------------------------

### load packages
library(foreign)
library(tidyr)
library(dplyr)
library(maptools)
library(gridExtra)
library(ggmap)
library(data.table)
library(raster)
library(rgdal)
library(sp)
library(grid)
library(gridExtra)
library(lmtest)
library(sandwich)
library(sf)
library(geosphere)

### empty environment
rm(list=ls())

### Set working directory
path <- "" # Define path here
setwd(path)



# Write function to compute fragmentation ---------------------------------

### write a function
fragmentation <- function(year){
  # load shapefile
  if(year<1300){
    shp <- readOGR(paste(path,"/Raw data/Fragmentation Abramson 2017 shapefiles/",
                         year,"_poly.shp",sep=""))
  }
  if(year>=1300){
    shp <- readOGR(paste(path,"/Raw data/Fragmentation Abramson 2017 shapefiles/poly_",
                         year,".shp",sep=""))    
  }
  
  # implement name changes (code from Abramson 2017)
  shp@data$Name<-ifelse(shp@data$Name =="Balearic islands" ,"Balearic Islands",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Albanians" ,"Albania",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Almavorid_4" ,"Almoravid_4",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Almafi","Amalfi",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Bresancon","Besancon",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Stratsbourg","Strasbourg",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Russians" ,"Russians_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Correggio" ,"Corregio" ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Wurrtemberg","Wurttemberg",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Ertfurt" ,"Erfurt" ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Regensberg" ,"Regensburg"  ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Ullster","Ulster",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="The  Holy Roman Empire","The Holy Roman Empire",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Fatimad_1","Fatimid_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Fatimad_2","Fatimid_2",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Fatimad_3","Fatimid_3",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Fatimad_4","Fatimid_4",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Fatimad_5","Fatimid_5",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Fatimad_6","Fatimid_6",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Kevian Rus","Kievan Rus",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Lithuania","Lithuanians",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Polaban Slavs","Polabian Slavs",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Poblabian Slavs","Polabian Slavs",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Pechengens","Pechenegs",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Kingdom of Burgandy","Kingdom of Burgundy",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Leistner","Leinster",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Noresmen","Norsemen",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="South_Slavs","South Slavs",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Strachclyde","Strathclyde",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaa_1_4_1","Ummayad_1_4_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_1_1","Ummayad_1_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_1_2","Ummayad_1_2",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_1_3","Ummayad_1_3",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_1_4","Ummayad_1_4",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_1_4_1","Ummayad_1_4_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_1_4_2","Ummayad_1_4_2",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_2","Ummayad_2",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_3","Ummayad_3",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_3_1","Ummayad_3_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_3_2","Ummayad_3_2",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_4","Ummayad_4",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_4_1","Ummayad_4_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_4_1_1","Ummayad_4_1_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Ummayaad_4_1_1","Ummayad_4_1_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_4_1_2","Ummayad_4_1_2",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_4_2","Ummayad_4_2",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_5","Ummayad_5",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaard_3_1","Ummayad_3_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayyad_3_2","Ummayad_3_2",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Geona" , "Genoa",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Almovorads", "Almoravids",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Balearic" | shp@data$Name ==         "Balearic islands"  , "Balearic Islands",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Hammadid", "Hammadids",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Norseman", "Norsemen",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Polostk", "Polotsk",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Munstert", "Munster",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Alhomad" | shp@data$Name ==  "Almohad" , "Almohads",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Bresica"  , "Brescia"  ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==  "Bulgarians"  ,  "Bulgaria"   ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Catile" | shp@data$Name == "Castille", "Castile" ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==    "Cremoa"   ,   "Cremona"    ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==    "Dauhphine"   ,    "Dauphine"  ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Felte and Belluno"  | shp@data$Name == "Feltre amd Belluno",  "Feltre and Belluno"  ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==   "Franche Compte" ,    "Franche Comte" ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==   "Kingdom of the Sword" ,     "Knights of the Sword"  ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==   "Kingdon of Naples"| shp@data$Name == "The Kingdom of Naples" ,  "Kingdom of Naples"  ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==   "Pistoa"   ,     "Pistoia"  ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==  "Portgual"     ,    "Portugal"   ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==   "The Republic of Novgorod"      ,   "Republic of Novgorod"   ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==    "Trevisto" ,   "Treviso"   ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Entense","Estense",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Freising","Friesberg",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="GIlan","Gilan",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Granada","Grenada",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Hapsburgs" ,"Habsburgs",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Il Khanid Mongols", "Il-Khanid Mongols",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Knights Hospitaller", "Knights Hopitaller",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==  "Luxembourg", "Luxembourgs",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==  "Mongol Empires", "Mongols",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==  "Mongol Empire", "Mongols",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==  "Vulga Bulgars", "Volga Bulgars"  ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==  "Wittelsbachs", "Wittlesbachs",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==  "Wurtemburg", "Wurttemberg",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==  "Hainaut", "Hainault",as.character(shp@data$Name))
  
  shp<-subset(shp,!is.na(shp@data$Name))
  shp@data$X<-as.numeric(as.character(shp@data$X))
  shp@data$Y<-as.numeric(as.character(shp@data$Y))
  shp@data$Area<-as.numeric(as.character(shp@data$Area))
  shp@data$logArea<-log(shp@data$Area)
  shp<-subset(shp,shp@data$Area > 6)
  
  # convert to sf file
  shp2 <- st_as_sf(shp)
  
  # adjust coordinate systems
  st_crs(shp2) <- "+proj=aea +lat_1=43 +lat_2=62 +lat_0=30 +lon_0=10 +x_0=0 +y_0=0 +ellps=intl +units=m +no_defs"
  shp_km <- st_transform(shp2, 32632)
  
  # determine centroids
  cent <- centroid(shp)
  cent <- SpatialPoints(cent, proj4string=CRS(as.character("+proj=aea +lat_1=43 +lat_2=62 +lat_0=30 +lon_0=10 +x_0=0 +y_0=0 +ellps=intl +units=m +no_defs")))
  
  # convert and create buffer
  cent <- st_as_sf(cent,coords=c("longitude.y","latitude.y"))
  cent_km <- st_transform(cent, 32632)
  cent_buffer100 <- st_buffer(cent_km, 100000)
  cent_buffer250 <- st_buffer(cent_km, 250000)

  # create output file
  out <- data.frame(matrix(nrow=dim(shp@data[1]), ncol=8))
  colnames(out) <- c("Name", "year", "frag100N", "frag100Meanarea", "frag100Medianarea",
                     "frag250N", "frag250Meanarea", "frag250Medianarea")
  out[,1] <- as.character(shp@data$Name)
  out[,2] <- year
  
  
  
  # intersection
  for(j in 1:nrow(out)){
    # define the intersections
    int100 <- st_intersection(cent_buffer100[j,], st_buffer(shp_km, 0))
    int250 <- st_intersection(cent_buffer250[j,], st_buffer(shp_km, 0))

    # compute N, mean area, median area
    out[j,3] <- dim(int100)[1]
    out[j,4] <- mean(int100$Area, na.rm=T)
    out[j,5] <- median(int100$Area, na.rm=T)
    out[j,6] <- dim(int250)[1]
    out[j,7] <- mean(int250$Area, na.rm=T)
    out[j,8] <- median(int250$Area, na.rm=T)  
  }
  
  # output
  return(out)
}



# Implement the function --------------------------------------------------

### Let's try this
# start with year = 1100
out1100 <- fragmentation(year = 1100)
fragmentationDATA <- out1100

# add all other years
for(y in seq(1105,1790, 5)){
  print(y)
  fragmentationDATA <- rbind(fragmentationDATA, fragmentation(year = y))
}

# drop repeat observations
fragmentationDATA2 <- subset(fragmentationDATA, !(Name == "Milan" & year == 1130))
fragmentationDATA2 <- subset(fragmentationDATA2, !(Name == "Milan" & year == 1135 & frag250N == 21))
fragmentationDATA2 <- subset(fragmentationDATA2, !(Name == "Milan" & year == 1170 & frag100N == 7))
fragmentationDATA2 <- subset(fragmentationDATA2, !(Name == "The Holy Roman Empire" & year == 1270 & frag100N == 9))

# save
write.dta(fragmentationDATA2, file="fragmentation.dta")


